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1 .  INTRODUCTION 

There  are  three  versions  of  the  finite  element  method.  The  classical  h- 

~  - '  'y  /  ^  v  — ' 

version,  which^  achieves/bh»  accuracy  by  refining  the  mesh  while  using  low  degrees 

p  of  the  elements,  usually  p  -  1,2.  The  p-version  keeps  the  mesh  fixed  and  the 

w 

accuracy  is  achieved  by  increasing  the  degree  p.  The  h-p  version  combines  both 


approaches.  jr _ .... 

The  p  and  h-p  versions  are  new  developments The  p-version  was  implemented  at 
Washington  University  in  St.  Louis  in  an  experimental  code  called  COMET-X  in  the 
middle  of  1970.  The  essential  part  of  the  code  were  hierarchic  elements.  This  type 
of  elements  was  first  considered  by  Zienklewicz,  Irons,  Scott  and  Campbell  [1970]  in 
conjunction  with  joining  finite  elements  of  different  polynomial  degrees.  Hierar¬ 
chic  C°  elements  were  then  described  by  various  authors,  e.g.  Peano  [1975],  Katz, 
Peano  and  Rossow  [1978],  Szabo  and  Peano  [1983],  Zienkiewicz,  Gago  and  Kelly  [1983]. 
The  cohesive  description  of  the  p-version  has  been  given  in  Szabo  [1979]. 

The  first  theoretical  analysis  of  the  p-version  was  given  in  Babuska,  Szabo  and 
Katz  [1981];  and  in  Babuska  and  Szabo  [1982].  The  performance  of  the  h-p  version 
was  first  theoretically  studied  in  Babuska  and  Dorr  [1981].  For  theoretical  analy¬ 
sis  of  the  p-version  in  3-dimensions,  we  refer  to  Dorr  [198H]  and  Dorr  [1986].  Ad¬ 
ditional  recent  results  are  mentioned  below. 

For  the  implementational  aspects  of  the  p-version,  we  refer  to  Szabo  [1985], 
Szabo  [1986]  and  Szabo  [1986a]. 

The  p  and  h-p  versions  for  two  dimensional  problems  were  implemented  in  the 
commercial  system  PROBE  by  Noetic  Tech.,  St.  Louis  with  first  release  in  1985,  and 
the  second  one  in  1986  (computations  in  the  present  paper  are  made  by  PROBE).  The 
three  dimensional  finite  element  code  FIESTA  having  some  p-version  capabilities  was 
developed  at  ISMES  (instituto  Sperimentali  Modelli  e  Strutture)  in  Bergamo,  Italy, 
and  since  early  1980  this  program  has  been  available  in  USA.  A  new  implementation 
for  three  dimensional  applications  on  Cray  computers  was  begun  by  the  Aeronautical 
Research  Institute  of  Sweden  (Glygtekniska  Forsoksanstalten  FFA) .  The  p  and  h-p  ver¬ 
sions  are  used  in  the  industry  today. 
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Although  the  p  and  h-p  versions  of  the  finite  element  method  are  relatively  new 
developments,  many  basic  results  are  available.  The  aim  of  -fe-ho  present  paper  is  to 
give  a  survey  of  the  basic  available  results  and  directions  for  further  development. 
The  paper  tries  to  survey  basic  theoretical  implementational  and  computational  as¬ 
pects  of  the  method  as  of  today.  / _ _ _ _ — - 


2.  FINITE  ELEMENT  METHOD  AND  THE  APPROXIMATION  PROBLEM 

Let  B(u,v)  be  a  bilinear  form  defined  on  H1  *  H2,  where  H1  and  H2  be 
reflexive  Banach  spaces  equipped  with  the  norm  |*|.|  and  |*|2»  respectively.  Let 
further  F  £  H2,  i.e.  F  be  a  linear  functional  on  H2. 

By  the  problem  (B,F}  we  denote  the  problem  to  find  uQ  £  H1  so  that 


B(u0lv)  -  F(v)  (2.1) 

holds  for  all  v  £  H2. 

If  the  bilinear  form  B(u,v)  is  continuous  and  satisfies  the  so  called  inf- 
sup  condition  (see  BabuSka,  Aziz  [1972],  Ch.  5),  then  the  problem  IB.F)  has  unique 
solution. 

Let  now  Sj  €  ^  €  H2.  Then  the  finite  element  problem  {B,F,S1tS2}  is 

to  find  the  finite  element  solution  uc  €  S,  such  that 

S1  1 

B(ug^ ,  v)  -  F(v)  -  B( Ug, v )  (2.2) 

holds  for  all  v  £  S2. 

If  the  bilinear  form  B(u,v)  satisfies  the  inf  sup  condition  on  S1  *  S2, 
then  exists,  is  unique  and 


where 


|ug  -  uQ| 1  $  C(S1,S2)Z(u0,H1,S1)  (2.3) 

Z(u0,H1,S1)  -  inf  |uQ  -  o)J  1 .  (2.H) 


For  detail,  see  BabuSka,  Aziz  [1972,  Ch.  6]  and  Arnold,  BabuSka,  Osborn  [1985].  We 
will  assume  that 


C(S1,S2)  $  D 


(2.5) 


where  D  is  independent  of  S1,S2  and  hence  the  norm  of  the  error  e  *  us  “  uo 
of  the  finite  element  solution  is  completely  given  by  Z(uq,H1,S1) 

Remark  2.1.  We  do  not  need  necessarily  that  (2.5)  holds.  Nevertheless,  as¬ 
sumption  (2.5)  simplifies  our  conclusions. 

Remark  2.2.  The  condition  (2.5)  is  satisfied,  for  example,  if  H,  •  H2, 
B(u,v)  -  B(v,u),  S1  -  S2  and 


This  condition  is  satisfied  for  selfadjoint  positive  definite  problems  as  in  the 
theory  of  elasticity,  etc. 


A 


The  exact  solution  uQ  is,  of  course,  not  known.  Nevertheless,  we  will  assume 
that  it  is  a-priori  known  that  uQ  €  K  c  Hj ,  where  K  is  a  certain  set  called  the 
solution  set  which  is  compact  in  H1 .  We  define 

Z^.HlS,)  -  sup  ZCu.Ht.S,)  (2.6) 

u  (K 

which  characterizes  the  error  under  the  assumption  that  we  know  only  that  the  solu¬ 
tion  u0  €  K. 

Remark  2.3.  A  typical  example  is  that  H1  «  and 

<  •  (u  |  |u|  .  *  1,  k  >  1). 

H  (n) 


This  choice  leads  then  to  the  classical  error  estimate  of  the  error  of  the  finite 
element  solution  (h-version): 


le  I  ,  <  Ch>  | 


See,  e.g.,  BabuSka,  Aziz  [1972,  Ch.  1J].  A 

There  are  many  results  available  concerning  the  characterization  of  Z(K, H-, ,  S1 ) , 
the  best  selection  of  S1  of  dimension  n,  etc.  For  an  excellent  abstract  treat¬ 
ment  and  survey  of  available  results,  we  refer  to  Pinkus  [1985]. 

The  space  S.|  in  (2.6)  is  fully  characterized  by  the  finite  element  method, 
its  h,  p  or  h-p  versions.  The  set  K  is  characterized  by  the  class  of  problems  to 
be  solved.  Hence,  the  performance  of  the  finite  element  method  relatively  to  the 
solution  set  K  is  described  by  Z(K,H,S)  which  will  be  in  the  center  of  our 
interest.'  Of  course,  others  aspects  are  also  essential  for  the  assessment  of  the 
performance  of  the  finite  element  method.  They  will  also  be  described  in  this 
survey. 


3.  THE  MODEL  PROBLEM 


The  performance  of  the  method  depends  strongly  on  the  class  of  problems  for 

which  it  is  used.  As  said  in  Section  2,  the  performance  is  directly  related  to  the 

solution  set  K  under  consideration.  We  will  concentrate  here  especially  on  the 

class  of  problems  which  are  characterized  by  the  piecewise  analytic  input  data. 

2 

Let  OCR1-  be  a  bounded  domain  and  its  boundary  be  a  piecewise  analytic 

n 

curve  T  -  U  ri  where  are  (closed)  arcs  with  the  end  points  AiPAi+1,  l  » 

i  «1 

1,...,n  (An+1  -  A-j ) .  An  example  and  the  notation  is  shown  in  Fig.  3.1. 


3h  ojn  ^rZA, 


Fig.  3.1-  The  scheme  of  a  domain  with  piecewise  analytic  boundary. 

By  Aj,  i  *  we  denote  the  vertices  of  ft  and  by  ,  i  »  1,2, ...n  the 

internal  angles.  We  will  not  exclude  the  case  when  the  internal  angle  uj  -  2n.  This 
case  is  very  important  in  practice  (cracks)  when  two  arcs  (fully  or  partially) 
coincide . 

Let  fD  =  U  and  rN  -  f  -  fD  be  the  Dirichlet  and  the  Neumann  boundary, 


respectively.  We  shall  be  interested  in  solving  the  problem 

-Au  +  u  »  f  on  ft 
u  -  h  on  fp. 


3n  “  8  on  V 


(3.1a) 

(3.1b) 


(3.10 


We  will  cast  the  problem  (3-D  (for  h  -  0)  into  the  form  of  a  {B, F}  problem.  To 
this  end  let  H  -  H-j  »  H2  -  Hp(ft)  -  {u  €  H1  ( ft)  |  u  -  0  on  rD)  where  by  H1  ( ft)  we 
denoted  the  usual  Sobolev  space  of  functions  with  the  square  integrable  first  deriv¬ 


atives.  Let 


F(v )  -  j  fv  dx  dy  ♦  J  gv  ds. 

ft  r„ 


(3.2a) 


(3.2b) 


If  h  4  0,  then  as  usual  we  write  u  «  z  +  w  with  w  £  H  ( ft) ,  w  =  h  on  rD  and 
z  €  Hp( ft)  being  the  solution  of  the  {B,F^}  problem  with  Fj  being  properly  ad¬ 
justed  F. 

The  model  problem  { B, F}  satisfies  the  conditions  listed  in  Section  2  provided 
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that  f,  h  and  g  satisfy  some  mild  conditions  as,  for  example,  f  €  L2(fl),  g  € 
L2(rN)  and  h  €  H1  (1^),  i  fc  V  and  h  continuous  on  fD- 

The  illustrative  numerical  computations  presented  in  the  next  sections  are  re¬ 
lated  to  the  two  dimensional  elasticity  problem,  i.e.,  for  the  (strongly  elliptic) 
system  of  two  partial  differential  equations  of  second  order  instead  of  the  simple 
model  problem  mentioned  above.  The  elasticity  problem  has  very  similar  property  as 
the  introduced  model  problem  but  has  larger  practical  importance. 

The  finite  element  solution  (for  h  ■  0)  is  characterized  by  the  selection  S1 
-  Sj>  -  S  c  hJj  and  all  conditions  inclusive  condition  (2.5)  in  Section  2  are  satis¬ 
fied.  If  h  4  0  and  h  is  not  a  trace  of  a  function  in  S,  then  we  replace  h 
by  h'  which  is  a  trace  of  a  function  in  S  and  consider  the  additional  error 
caused  by  this  replacing. 

Characterization  of  the  solution.  Set  K  relates  to  the  available  information 
about  input  data,  i.e.  the  information  about  T,  f,  g,  h.  We  will  assume  that 
r  is  piecewise  analytic,  f  is  analytic  on  Q,  g,  h  are  analytic  on  f  This 
assumption,  namely,  that  the  data  are  piecewise  analytic,  is  practically  always 
satisfied  in  the  problems  of  structural  mechanics. 

Remark  3-1.  In  our  illustrative  computations  we  will  also  present  the  results 
which  are  outside  of  the  above  mentioned  frame,  namely  that  g  is  a  Dirac  function 
(concentrated  load).  Such  an  example  is  well  taylored  for  our  illustrations,  but 
needs  more  refined  theoretical  analysis  which  will  not  be  addressed  here.  & 

Although  our  main  emphasis  will  be  on  the  problem  with  piecewise  analytic  input 
data,  we  will  also  mention  the  results  for  the  more  usual  solution  set  K  as,  for 

example,  K  -  {u  |  |u|  SI),  etc. 

HK(£2) 

Usual  assumptions  in  the  regularity  theory  of  the  differential  equations  of  el¬ 
liptic  type  are  based  on  the  theory  of  Sobolev  spaces  of  finite  order,  i.e.  f  € 

H  (Q),  g  €  H  (D,  etc.,  and  often  the  boundary  of  the  domain  is  assumed  to  be 
smooth  (i.e.  not  piecewise  smooth).  Such  assumptions  are  not  sufficiently  realistic 
in  applications.  Either  they  are  too  restrictive  (smooth  domain  G)  or  too  general 
(f  €  Hk(G))  _ 

To  further  simplify  our  exposition,  we  will  assume  that  fl  is  a  polygon.  We 
will  make  some  remarks  about  the  general  case. 

Remark  3.2.  We  will  also,  as  illustrative  example,  deal  with  one  dimensional 
analog  of  our  model  problem,  namely  the  problem 

-u"  -  f,  x  €  (0,1)  -  I  (3-2a) 

u(0)  -  u(7)  - 

with  f  such  that  the  exact  solution  Uq(x)  is 


0 


(3.2b) 


where 


(x  -  O 


x  -  e 

o 


and  a,  b  are  such  that  (3.2b)  is  satisfied. 


for  x  >  C 

for  x  <  5 

Obviously  uQ  F  H1 ( I) . 


A 


4.  CHARACTERIZATION  OF  THE  SOLUTION  SET 

As  it  was  said  earlier,  the  solution  set  K  describes  the  solutions  of  the 
class  of  problem  to  be  solved.  The  performance  of  the  method  is  then  directly  re¬ 
lated  to  this  set. 


Let  B_  =  (61 . Bn)  be  a  n-tuple  of  real  numbers  0  <  8^  <  1 ,  1  £  i  £  n.  For 

any  integer  k  >  0  we  shall  write  8  +  k  «  ( B-| +k ,  B2+k , . . . ,  Bn  +  k).  By  r^x),  i.  = 
we  shall  denote  the  Euclidean  distance  between  x  €  £2  and  the  points  Bi 

n  e. 


Q. 

n 


,n.  We  denote  then 


i  ■  1 . 

B+k 

nv  (x).  The  points 

ierti 


®8+k(x) 


TT  »*!*(*) 

i-1 


and  $ 


B+k 


(x) 


ces 


or  outside. 


Bi  could  be  located  at  the  boundary  of  n,  e.g.  in  the 
They  also  can  be  absent,  but  we  will  not  elaborate  on 


this  case  in  this  paper. 
Define  now 


^  -  (u  €  Hp(fl)  |  1/  | D°u| 

a 


2  *2 


B+k -2 


(x)dx  dy Y*  < 


Cdk  2 (k-2) ! , 


k  -  2, 3, . . . , | a|  -  k.  d  >  1,  d  independent  of  k}. 

As  usual,  we  denoted  a  -  (aj,a2).  |a|  -  a1  +  a2,  ai  £  0,  i  »  1,2,  integers  and 

d“u  _  alalu 

a1  a2  * 
dx'dy 

n 

The  functions  belonging  to  are  analytic  on  55  -  £  Bj.  If  Bi  €  3Q,  then 

i=l 

they  have  singular  behavior  in  the  neighborhood  of  ,  and  the  character  of  the 
singularity  is  given  by  Bi  and  d. 

It  has  been  shown  in  Babuska,  Guo  [1986]  that  if  the  domain  Q  is  a  polygon, 

»  Aj^  (i.e.  Bi  are  the  vertices)  and  functions  f,  g,  h  are  analytic  on  55 
and  rif  respectively,  then  the  solution  of  the  problem  (3-1)  belongs  to  for  a 

properly  chosen  constants  B,  C,  d.  The  case  B^  ^  55  characterizes  the  solu¬ 
tions  with  the  singularities  outside  of  Q,  e.g.,  when  the  domain  has  circular  arcs 
and  h  -  0.  This  case  describes  also  well  the  case  when  the  natural  domain  of  the 
analyticity  of  the  solution  contains  55. 

The  set  K,  obviously  belongs  to  the  family  of  countably  normed  spaces.  For 
more  about  this  family,  we  refer  to  Gelfand,  Shilov  [1964]. 


Let  us  now  introduce  the  more  standard  family  of  solutions  sets 


K2  =  {u  €  h!(Q)  |  |u|  <  C,  k  >  1 } 

K 

h  !(n) 

<3  =  [uf  h’CO)  |  |u|  k  <  C,  k2  > 

H  2(fl) 

=  {u  €  Hp(fl)  |  u  -  r1<Xi|log  ^  |  1  ^  (6i  )x*  (rt ) } 

where  (r^.0^)  are  the  polar  coordinates  with  the  origin  in  the  vertex  Ai ,  ai  > 

0  noninteger,  'P i ( 6^ )  is  a  c“  cut  off  function. 

The  motivation  of  the  solution  set  Kj ,  j  =  2,3,4  is  that  the  solution  u  of 
(3-1)  can  be  written  in  the  form 

U  =  U1  ♦  U2  +  Uj  ( 4. 1 ) 

where 

u-j  €  K2,  U2  ^ 

and 

u3  ■  l  ci,juiJ]*  uiJ]  C  K4* 


Functions  u1  and  u^  satisfy  the  homogenous  Dirichlet  (essential)  conditions, 
while  u2  relates  to  the  nonhomogeneous  Dirichlet  conditions.  The  restriction  k2 
>  3/2  has  been  made  for  simplicity  only  and  can  be  replaced  by  k^  >  1 .  For  the 
theory  leading  to  the  form  (4.1)  we  refer  to  Kondrat'ev  [1967],  Kondrat'ev,  Olejnik 
<1983]  and  Grisvard  [1985]. 


Remark  4,1.  We  restricted  ourself  to  the  problem  (3.1)  only.  The  practically 
important  case  of  nonhomogeneous  materials  which  is  described  by  the  equation 


!_  a  +  i_  a  iH. 

3x  3x  3y  3y 


f 


with  a  being  a  piecewise  constant  on  the  domains  bounded  by  the  piecewise  analytic 
curves  can  be  handled  in  a  similar  way.  Analogous  situation  occurs  also  when  deal¬ 
ing  with  the  problem  of  elasticity.  & 


5.  THE  FINITE  ELEMENT  SPACE  S 


We  introduce  now  the  finite  element  spaces  we  will  deal  with  later. 

For  reasons  of  the  simplicity  of  the  exposition,  we  will  restrict  our  choices, 
but  our  numerical  example  will  also  present  more  general  cases. 

Let  M  -  {T}  be  a  family  of  meshes  T  -  in}  where  c  n  is  an  open  trian¬ 
gle,  called  element.  Let  hT  -  diam  xif  h(x)  -  max^  hT  ,  and  let  pT  be  the 


diameter  of  the  largest  ball  contained  in  x^.  We  will  assume  that  M  is  such  that 
for  any  xj  (  T  (  H 


i.e.  that  all  triangles  satisfy  the  minimal  angle  condition.  Further  we  shall  as¬ 
sume  that  2  =  U  Tj  and  that  any  pair  x,  ,  x,  €  T,  i  ^  j  has  either  an  entire 
x^T  j 

side  or  a  vertex  in  common,  or  their  intersection  is  empty.  The  number  of  elements 
of  t  will  be  called  cardinality  of  T  and  denoted  by  M(T) . 

Remark  5.1.  We  restricted  ourself  to  triangles  only.  The  results  we  will  pre¬ 
sent  are  valid  more  general,  e.g.  for  rectangles,  curvilinear,  triangles  and  rect¬ 
angles.  A 

In  addition  to  a  general  mesh,  we  will  introduce  two  special  families  of 
meshes ,  the  quasi uni form  and  the  geometric  mesh. 

a)  The  family  m7,,  1  £  Y  <  «,  of  quasiuniform  meshes:  There  exists  con¬ 


stant  y  such  that 


h(x)  »  max  h 


<  y  min  h 


holds  for  any  t  €  My. 

b)  The  family  Mg*k ( B1 , . . . , Bn) ,  Bj  $  2,  j  -  1 . .  0<q<1,  1  £  k  < 


®  of  geometric  meshes: 


Let  ti  €  Mg*k.  fl  l  B  -  0 

j  =  l  3 


m  0,  then  there  exists 


min  d(B. ,x. ) 
1<j<n  J 


k-1  l^q  <(i) 
q  -  h(xi) 


<  l^k 
q 


where  by  dCBj.x^)  we  denoted  the  euclidean  distance  between  Bj  and  x^ 

Obviously,  if  x  €  My',  then  M(T)  -  h“2(T)  and  the  cardinality  M(x)  can  be 

arbitrarily  large.  If  Bj  {  2,  j  -  1 . n,  then  for  any  we  have  a 

M(T)  <;  Q  where  Q  depends  on  2,  B,  q,  k.  If  Bj  €  32,  then  the  cardinality 
can  be  arbitrarily  large. 

Finally  we  denote  S(T,p)  -  {u  €  Hp( £2)  |  u|  is  a  polynomial  of  degree  £ 
p  for  any  xi  €  T}.  By  N(S,T,p)  -  dim  S(T,p)  we  denote  the  number  of  degrees  of 
freedom . 

If  T  €  Mj,  then  N(S,  ,p)  -  p2h“2(  ).  If  x  €  g’k,  Bi  f  a,  then 

N(S,T,p)  »  2  and  If  T  €  M^,k  with  B^  €  32,  then  N(S,  ,p)  «  p2|  log^ n(T)  |  , 
Remark  5.2.  We  assumed  only  triangular  elements.  For  the  rectangular 
elements,  instead  of  polynomials  of  the  total  degree  p,  we  use  polynomials  of 
degree  p  in  every  variable.  In  the  case  of  the  curvilinear  elements  we  use  the 
standard  mapped  polynomials.  A 


Remark  5.3.  We  mentioned  only  the  case  when  the  degree  p  is  uniform,  i.e.  is 
the  same  for  all  elements.  The  available  theory  covers  also  the  nonuniform, 
selective  choice  of  the  degrees  p.  A 

Remark  5.1*.  The  meaning  of  the  finite  element  subspaces  is  especially  clear  in 
the  one  dimensional  setting.  For  example,  in  the  case  of  the  geometric  mesh  with  A 
-  0,  then  nodal  xt  points  are 

Xi  =  qM_i,  i  =  0,1 . M 

where  q  <  1  is  called  the  ratio  of  the  mesh.  A 


6:  THE  BASIC  APPROXIMATION  THEOREMS 

We  will  mention  here  some  basic  approximation  for  H1  =  H^(ft),  i.e.  we  will 
study  the  magnitude  of  Z(/Cj , H 7 ( Q) ,  S(T,p)).  Numerical  illustration  will  be  given 

in  Section  7. 

First  let  us  consider  any  fixed  mesh  Tq  and  consider  the  p-version. 


Theorem  1 . 

Z(K2,Hi(0),S(T0,P))  <  Cp-'(kr1) 


(6.1 ) 


holds  for  any  mesh  TQ  (satisfying  (5.0).  For  the  proof,  see  Babuska,  Suri  [1985]. 
Theorem  2. 

Z(Kjj,  Hp(  fl) ,  S(7"0,  p ) )  <;  C|  log  p|V2a  -  C|  log  N|  Va.  (6.2) 


For  the  proof,  see  Babuska,  Suri  [1985]. 

Theorem  3.  t  K  -  and  Bj  ^  Q,  i  »  1,...,n.  Then 

»/ 

Z(f(1,Hj)(n),S(T0,p))  1  Ce'ap  -  Ce'aN*  (6.3) 


where  a  depends  on  B^ ,  i  »  1,...,n  but  is  independent  of  p.  For  the  proof, 
see  Guo,  BabuSka  [1986]. 

So  far  we  assumed  that  the  essential  boundary  conditions  are  homogeneous.  Let 
us  discuss  now  the  case  K ^  when  the  Dirichlet  conditions  are  not  homogeneous.  Let 
r  c  rD  be  any  side  of  an  element  x  (  T  with  the  end  points  P^,  i  =1,2  and 
let  <p  be  defined  on  T  so  that 

a)  (p  is  a  polynomial  of  degree  p 

k  p 

b)  cpCPj)  =  uCPp,  i  =  1,2,  u  €  K3  (because  u  (  Kj,  u  €  H  *(Q),  k2  > 

3/2,  u(Pj)  is  well  defined) 

c)  f  ds  =  j  u'ip'  ds  for  any  ip  being  polynomial  of  degree  p  and 

F  f 

<P(P.)  =0,  i  =  1,2. 


Now  we  replace  the  boundary  condition  on  fD  by  <p,  i.e.  us  =  ip  and  then  we 


Theorem  4. 


Z(fC3,H1(n),S(x0,p))  <  Cp 


(k2-1) 


For  the  proof,  see  Babuska,  Suri  [1985]. 

Remark  6.1.  In  the  one  dimensional  case  we  can  prove  more  exact  theorems.  Let 
us  consider  the  problem  (3.2)  with  the  solution  (3-3)  and  the  case  of  one  element 
only.  Then  we  have 


Theorem  5.  Let  Hg(I)  =  {u1  €  L2(I),  u(0)  =  u(1)  =  0}  and  | u | _  1 


|u'Il2(I)-  Then 


H0(D 


a)  if  5=0  then 


Z(uq,  Hq(  I) ,  S(p) )  =  CQ(a)  j-  (1+0(1))  ,  [p  -  »)  (6.5a) 


Cn(o) 


ar(g)  |sin  ira| 

it/2a-1 


b)  if  i|>  <  0,  then 


where 


Z(uq,  Hg(  I) ,  S(p) )  -  Cl(a)[^-)a  1  ^  iTOl-^)) 


o  >  0.  r-  {EJLlJEL  t  c  (a)  _  ar(a)  [sin  ira| 

✓TT  *  1  /?  2a“4 


(6.5b) 


(6.6a) 


(6.6b) 


the  form  0(— )  is  uniform  with  respect  to  (  (  -e,  e  >  0  and  we  have  the  estima- 
P  2  1 

tlon:  if  0  <  £  1  -  —  ,  then 


Z( u0, Hq( I) , S(p ) ) 


if  1  -  —  r2  <  1 ,  then 
P 


1  rP  1  -a  r  1  ,,  2,°-^ 

/1-r2  P  p 


p+1-a  .  a-l4 

r  1  .  , .  *  i 


Z(u0.Hj(I).S(p))  »  r  j  [-jUr  »  (1-r)  ] 


a  ~4  1  a -'4 

p  2  p  2 


(6.6c) 


(6.6d) 


where  equivalency  constants  depend  only  on  a. 


c)  if  0  <  5  <  1,  then  there  exists  a  constant  C  >  0  depending  only  on  a 


such  that 


where 


g(h,p,Y)  =  max(|log  h|\|log  p|Y). 

For  the  proof,  see  Babuska,  Suri  [1968a]. 

Finally,  let  us  address  the  case  of  the  h-p  version  on  the  geometric  mesh. 
Theorem  8.  Let  K  -  and  ler  B^  be  vertices  of  Q.  Then  there  exists  a 
mesh  and  the  degree  p  (dependent  on  the  mesh)  such  that 

3— 

Z(^  ,hJ(0),S(t,p))  <  Ce_ct™  (6.9) 

where  a  is  independent  of  N. 

For  the  proof,  see  Guo,  Babuska  [1986]. 

Remark  6.2.  If  the  degrees  of  the  elements  are  uniform,  then  the  class 
can  be  extended.  See,  Guo,  Babuska  [1986c].  A 

Remark  6.3.  In  one  dimensional  setting  much  more  can  be  said.  We  will  discuss 
the  case  5=0. 

Let  us  first  ask  the  question  about  the  lower  bound  of  the  error  Z(u0,Hq(I), 
S(T,£))  among  all  meshes  and  all  degree  distributions.  The  answer  is  given  in 

Theorem  9.  _ 

Z(uQ, H^j(  I) , S( T »£_) )  *  C(a)  —  -■  qg(a~1/;)N  (6.10) 

where 

90  -  1^-1  )2. 

For  the  proof,  see  Gui ,  Babuska  [1986]. 

Let  us  now  consider  the  geometric  mesh  with  the  ratio  q  and  various  degree 
distribution.  We  have 

Theorem  10.  As  N  ♦  «  the  optimal  degree  distribution  tends  to  be  linear  with 
the  slope 


3n  -  (a-1/,) 


vi  ln  q 


2'  In  r  ’ 


(6.11a) 


ft 


(6.11b) 


I 


This  means  precisely  that  for  each  i  =  1,2,... 

(M  )  (M  ) 

lim^PMN-i  "  PMN-i-1 ^  ”  30* 

For  the  proof,  see  Gui,  Babuska  [1986]. 


For  the  error  estimates  with  linear  degree  vector  of  slopes  S  we  have  the  follow- 


Theorem  1 1 .  For  the  geometric  mesh  with  the  ratio  q  combined  with  a  linear 


degree  vector  of  the  slopes  s  we  have 
a)  if  s  >  sn,  then 


Z(u0,Hq(I)  ,S(T,p_))  »  C(a,q,s)q 


(a-‘/,)/2N/s 


(6.12a) 


b)  if  s  <  sn,  then 


Z(u0, Hq( I) , S(T,p_) )  -  C(a,q,s)rV 


c)  if  s  =  Sr.,  then 


Z(u0,Hq(I),S(T,p))  -  C(a,q)e 


-/( a-l/2) N  /21n  q  In  r 


(6.12b) 


(6.12c) 


where  r  =  1  and  sn  =  (a-l4)  t— — ■ —  is  the  optimal  slope  in  the  sense  that  the 
1=^  °  ln  r 

exponential  rate  attends  maximum.  Furthermore,  the  optimal  geometric  mesh  and  linear 


degree  vector  combination  is  given  by 


(/2  -  1)‘ 


2a  -1. 


In  this  case 


Z(u0,Hj(I),S(T,p))  -  C(a)[(/2  -l)2]^®  . 


(6.13) 


(6.14) 


In  (6.12)  -  (6.16)  the  equivalence  constants  depend  on  (a,q,s),  and  a,  respec¬ 
tively. 

For  the  proof,  see  Gui ,  Babuska  [1986]. 

Remark  6.5.  For  the  optimal  combination  of  geometric  mesh  and  linear  degree 
vector,  the  estimate  can  be  written  as 


,,  ul  t  .I  ..  -1 .762/(a-l/JN 

2(uqi Hq ( I) , S( t  ,p ) )  *  e 


(6.15) 


and  we  have  seen  in  Theorem  9  that  this  exponential  rate  of  convergence  is  the  best 
possible  one.  t 

Remark  6.5.  We  discussed- in  Remark  6.2  the  case  of  linear  degree  vector  of 
slopes  3.  Let  us  discuss  now  the  case  of  the  uniform  vector  of  degrees. 

Theorem  12.  For  the  geometric  mesh  with  the  ratio  q  combined  with  uniformly 
distributed  degree  p,  the  relation  between  the  optimal  choice  of  p  and  the  num¬ 
ber  of  element  M  in  the  mesh  is  asymptotically  linear,  i.e. 

p  »  SqM  (as  M  ♦  ®) 

with  Sq  being  the  same  as  in  Theorem  10. 

For  the  proof,  see  Gui,  Babuska  [1986]. 

The  analog  of  Theorem  11  is: 


§1 


8 


Theorem  1 3 •  For  the  geometric  mesh  with  ratio  q  and  the  uniformly  distrib¬ 
uted  degree  p  related  with  the  number  of  elements  M  by  p  =  sM,  we  have 
a)  if  s  >  Sq,  then 

(a-‘4)/N77 

Z(u0,  Hq(  I) , S(  T,p ) )  .  C(a.q.S)  ^ j—  (6.16a) 


(/sN)' 


b)  if  s  <  sn,  then 


Z(u0,Ho(I).S(T,p))  -  C(a,q) 


c)  if  s  =  s0,  then  one  gets  optimal  rate  of  convergence 


(6. 1 6b ) 


§ 

y 


where 


Z(u0, Hq{  I) ,  S(  T,p ) )  -  C(ct,q)  - 


-/ ( a-*/j) N  /In  r  In  q 


(a-‘4)ln  q 
In  r  ’ 


1  -Vq 

1+q  * 


min(2a-1 ,a) , 


The  optimal  combination  is  also  given  by 


(/2  -  1)‘ 


(6.16c) 


(6.1 6d ) 


(6.1 6e ) 


2a  -  1 


and  for  optimal  combination  we  get 


/(a-1/)N/2 


Z(uQ. Hg( I) , S(T,p) )  -  C(o) 


For  the  proof,  see  Gui ,  Babuska  [1986]. 


(6. 1 6f ) 


Remark  6.6.  We  can  also  interpret  the  optimal  h-p  version  with  uniform  p  as 
the  envelope  of  the  h  version  with  fixed  p.  In  this  case  the  for  N  -*■  »  the  mesh 


tends  to  be  geometric  with  the  ratio  q  «■  e_‘,/B  -  0.5820  and  the  relation  between 
the  degree  p  and  the  number  of  elements  M  tends  to  be  linear  with  p  -  ( 4/e2)  (a-*4)M 
-  0.541 3(a_14)M*  For  more  detailed  analysis,  we  refer  once  more  to  Gui,  Babuska 
[1986].  i 

Remark  6.7.  In  the  case  K  -  K i  Theorem  8  holds  when  we  restrict  ourself  to 
the  nonuniform  distribution  of  degrees.  The  result  holds  also  if  curvilinear  ele¬ 
ments  which  have  to  satisfy  certain  conditions. 

In  the  case  when  only  uniform  distribution  of  degrees  is  used,  then  these  con¬ 
ditions  are  weaker  than  in  the  general  case.  See  Guo,  Babuska  [1986c].  L 

Remark  6.8.  We  discussed  only  the  model  problem  (3-D-  The  results  hold  for 
more  general  equations.  For  the  higher  order  equations,  see  Guo,  Babuska  [1986b], 
and  Babuska,  Suri  [1986c].  In  the  case  of  systems  of  second  order  or  higher  order 
differential  equations,  there  is  much  broader  scope  of  essential  boundary  condi¬ 
tions.  For  their  treatment  we  refer  to  Babuska,  Suri  [1986b].  L 


Let  us  now  make  some  comments  to  the  theorems  we  mentioned  above. 

a)  Comparing  the  performance  of  the  finite  element  method  with  quasiuniform 
mesh  with  respect  to  the  number  of  degrees  of  freedom,  then  the  p-version  (with  few 
elements)  perform  better  than  the  h-version.  In  the  case  that  the  solution  has  sin¬ 
gularity  of  the  type  occurring  in  the  corner  of  the  domain,  the  rate  of  the  p-version 
is  twice  that  of  the  h-version.  For  the  h-version  there  is  a  classical  theorem  men¬ 
tioned  in  the  basic  books  about  finite  element  method 

inf  |u-u|  .  <  C(p)hn|u|  (6.17a) 

ui€S(t,p)  H  (Q)  H  (fl) 

n  =  min(p.k-l)  (6.17b) 

where  C(p)  is  not  specified  (more  precisely,  the  proof  indicate  that  C(p)  +  •  as 
p  •*  ®) .  This  leads  sometimes  to  the  (false)  statement  that  (6.17b)  indicate  that 
for  singular  solution  it  is  improper  to  use  higher  order  elements. 

b)  The  h-p  version  leads  to  the  exponential  rate  of  convergence  when  the  input 

2 

data  are  piecewise  analytic.  In  the  one  dimensional  case,  the  ratio  q  -  (/2  -1) 
is  the  optimal  one  Independently  of  the  strength  of  the  singularity.  In  the  two  di¬ 
mensional  case,  the  ratio  of  the  same  magnitude  seems  to  be  optimal  although  de 
tailed  theoretical  analysis  is  not  available  yet.  Some  numerical  evidence  will  be 
presented  below.  For  practical  reasons  q  -  0.15  is  recommended. 

c)  Although  we  mentioned  only  the  simple  model  problem,  the  results  hold  much 
more  generally.  We  mention  especially  the  elasticity  problem. 

d)  Very  important  problems  arise  in  relation  with  the  "locking"  problem  as  in 
the  elasticity  problems  with  near  incompressibility  (Poisson  ratio  v  -  %) .  It  has 
been  shown  by  Vogelius  [1983]  that  the  p  (and  h-p  version)  is  not  influenced  by  the 
locking  problem  and  solves  reliably  the  elasticity  problems  with  nearly  incompress¬ 
ible  materials  without  any  difficulties.  See  also  Babu3ka,  Szabo  [1982]. 

7.  NUMERICAL  ILLUSTRATIONS 

In  thi3  section  we  will  present  numerical  illustrations  related  to  the  theorems 
we  mentioned  in  the  previous  section. 

Example  7.1.  Let  us  consider  the  plane  strain  elasticity  problem  when  fl  is 
the  L-shaped  domain  shown  in  Fig.  7.1. 

Let  us  assume  that  on  3Q  tractions  are  prescribed,  i.e.  rD  =  9.  We  will  assume 
that  the  solution  of  this  problem  is  the  displacement  vector  u  -  (u1 ,u2)  where 

u1  -  ijr  ra  [(<-Q(a+1))cos  a9  -  a  cos(a-2)9] 

(7.1 ) 

u 2  “  r°[( <+Q(u+1 ) )sin  a9  -  a  cos(a-2)9] 

where 


Fig.  7.1.  L,-shaped  domain. 


a  -  0.544483737  Q  -  0.543075579. 

G  is  the  modulus  of  rigidity  and  <  -  3-4v,  where  v  -  .3.  The  sides  0A  and  0E 

are  traction  free.  The  solution  has  a  typical  singularity  at  0  and  is  the  first 

mode  of  the  stress  intensity  factor  solution.  Instead  of  the  norm  |*|  .  we 

H  (fi) 

will  be  interested  in  the  energy  norm  |»|  which  is  equivalent  to  the  norm 
1*1  ,  .  We  then  define  the  relative  error  lei-  »  |e |_/ |u I _.  . 

h\q)  e,r  e  e 

First,  we  will  consider  the  case  of  the  uniform  mesh  with  square  elements  shown 
in  Fig.  7.2. 

The  solution  u  can  be  interpreted  as  a  member  of  the  set  ,  with  B-|  »  0, 

K2  with  k2  -  a  -  e,  e  >  0  arbitrary  or  Kjj  with  Y  -  0.  Interpreting  the  so¬ 
lution  u  as  the  number  of  the  solution  set  Kj,,  the  estimate  (6.2)  gives 

hmin(a,p-a) 

u-us|  <  C  min[ha,  - ^ - 

P 

where  C  depends  on  a  but  is  independent  of  h  and  p. 

Fig.  7.3  shows  the  relative  error  lei.,  D  (for  different  p)  in  dependence 
on  h  (in  log|e|g  Rx|log  h|  scale).  We  also  show  the  slope  ha  in  the  figure. 

We  see  that  with  respect  to  h,  the  error  is  in  the  asymptotic  range  also  for 
moderate  p  and  h.  Figure  7.4  shows  the  error  in  dependence  on  p  and  differ¬ 
ent  h.  Because  of  the  size  of  the  computations  the  error  is  given  for  p  >  4  only 
for  h  -  \  ,  (for  p  -  4  and  h  -  1/10  the  number  of  degrees  of  freedom  N  -  5119). 

The  slope  p-2a  (2a  -  1.088)  is  apparent  only  for  p  >  3.  Figure  7.5  shows  the 
error  in  dependence  on  the  number  of  degrees  of  freedom  N  for  various  p.  Also  the 
performance  of  the  p-version  for  h  -  %  is  shown  in  Fig.  7.5.  We  see  that  the  p- 
version  is  more  effective  than  the  h-verston  and  that  the  theoretical  asymptotic 
slope  (shown  in  the  figure)  is  achieved  for  moderate  accuracy  and  N. 


Fig.  7.2.  The  scheme  of  the  uniform  mesh. 


Fig.  7-3.  The  relative  error  measured  in  the  energy  norm  in  dependence  on  h. 


Remark  7.2.  Let  us  illustrate  the  estimate  given  in  Theorem  5.  Introduce  the 


numerical  constant 


*P 


RA  - 

P 


where  is  the  right  hand  side  of  (6.5)  and  Ep  -  |e|-1 


H0(D 


Table  7.1  shows  the  error  Ep  and  the  numerical  constant  Rp  for  a  -  0.7, 

3A 


3.5  and  5-0.  We  see  that  for  a  smaller  R£  ♦  1  quicker  than  for  a  larger 
(i.e.  smoother  function). 


For  more  numerical  results  we  refer  to  Gui ,  Babuska  [1986], 


£ 

So  tar  we  addressed  the  performance  of  the  uniform  mesh.  Let  us  discuss  now 
the  performance  of  the  finite  method  on  the  geometric  mesh.  Fig.  7.6  shows  the  geo¬ 
metric  mesh  with  n  =  2  layers  and  the  ratio  .15. 


Fig.  7.4.  The  relative  error  measured  in  the  energy  norm  in  dependence  on  p. 


Fig.  7.5.  The  relative  error  measured  in  the  energy  norm  in  dependence  on  v. 

Fig.  7.7  shows  the  performance  of  the  p-version  for  various  number  n  of  layers 
in  log|e|g  R  *  log  N  scale.  Fig.  7.7  shows  also  for  comparison  the  performance  of 
the  method  on  the  uniform  meshes  shown  in  Fig.  7.2.  We  see  that  for  every  number  of 
layers  the  rate  is  n“-5441  where  N  is  sufficiently  large  (dependent  on  n),  and 
that  the  error  has  the  reverse  S  slope  form.  In  the  first  phase  the  accuracy  is 
exponential  (curved  down)  while  in  the  second  phase  the  accuracy  is  an  algebraic 
one.  Further  we  see  that  the  optimal  number  of  layers  and  the  degree  p  depends  on 


9 


the  required  accuracy.  We  will  return  to  this  question  later  in  the  discussion  of 
the  expert  system  for  the  h,p  version.  See  Section  11.  Fig.  7.8  shows  the  relative 
error  in  the  log|e|_  „  *  N1 ^  scales.  We  see  that  the  envelope  of  these  curves 
which  characterize  the  performance  of  the  h-p  version  is  nearly  a  straight  line  as 
it  has  been  theoretically  indicated  by  Theorem  8. 


TABLE  7.1.  The  error  and  the  values  of  the  numerical  constant 
in  the  one  dimensional  case. 


Sp 

nA  1 

% 

*5 

i|.  7^3  E-1 

0.9877 

1 .021 

0.2032 

3.627  E-1 

0.9967 

3.402  E-1 

4.335 

3.090  E-1 

0.9985 

3.093  E-2 

4.488 

2.756  E-1 

0.9992 

2.379  E-3 

1.940 

2.522  E-1 

0.9995 

4.760  E-4 

1.480 

2.344  E-1 

0.9996 

1.400  E-4 

1.300 

2.204  E-1 

0.9997 

5.154  E-5 

1 .208 

2.090  E-1 

0.9998 

2.210  E-5 

1.153 

1.994  E-1 

0.9998 

1.057  E-5 

1.118 

1.912  E-1 

0.9999 

5.495  E-6 

1 .094 

1.840  E-1 

0.9999 

3.053  E-6 

1.077 

1.777  E-1 

0.9999 

1.790  E-6 

1 .064 

1.722  E-1 

1.000 

8.999  E-7 

1.054 

1.671  E-1 

1.000 

6.978  E-7 

1.046 

1.626  E-1 

1.000 

4.585  E-7 

1 .040 

0  Q0225 


Fig.  7.6.  The  geometric  mesh  with  n  -  2  layers 
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Fig.  7.7.  The  performance  of  the  p-version  on  various  meshes  in  the  scale 
log|e|  „  *  log  N. 


The  envelope  together  with  the  optimal  pair  (n,p)  is  shown  in  Fig.  7.9, 
together  with  the  approximate  straight  line. 

The  p  and  h-p  version  has  various  properties  which  are  essentially  different 
from  the  h-version.  We  already  mentioned  the  robustness  of  the  p-version  with  re¬ 
spect  to  the  Poisson  ratio  (the  nearly  incompressible  material).  Let  us  discuss  now 
briefly  the  "pollution"  problem  which  describes  the  influence  of  the  locally  un¬ 
smooth  solution.  For  the  h-version  the  well  known  interior  estimate  shows  that  the 
influence  of  the  local  unsmoothness  of  the  solution  is  also  local.  For  the  p-version 
the  effect  is  not  so  local  as  for  the  h-version.  To  illustrate  this  effect,  let  us 
solve  the  plane  strain  elasticity  problem  <v  -  *3)  loaded  by  the  concentrated  load 
(Dirac  function)  as  shown  in  Fig.  7.10.  On  the  sides  BCDC'B'  the  tractions  are 
prescribed  so  that  the  exact  solution  is  the  classical  Businesque  solution  on  the 
half  plane  (hence  B'AB  side  is  traction  free).  The  energy  of  the  solution  is  infin¬ 
ite  and  hence  the  usual  theory  has  to  be  generalized.  We  will  be  interested  (be¬ 
cause  of  the  obvious  symmetry)  in  the  finite  element  solution  and  its  accuracy  on 
Q  (shadowed  in  Fig.  7.10).  As  before,  we  will  measure  the  accuracy  in  the  energy 
norm.  The  meshes  used  for  the  computation  are  shown  in  Fig.  7.11  (only  half  of  the 
domain  is  shown).  The  domain  S3  is  once  more  shadowed.  The  minimal  (pollution 
free)  error  on  0  is  the  error  of  the  finite  element  method  when  the  exact  trac¬ 
tions  are  prescribed  on  30. 

Fig.  7.12  shows  the  error  Ie I e( n)  for  the  meshe3  I_IV  and  the  "pollution 
free"  error.  We  see  that  the  p-version  for  any  shown  mesh  does  not  converge  al¬ 
though  for  three  layers  the  error  is  in  the  range  of  engineering  accuracy.  For 
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Fig.  7.11.  The  meshes  used  for  the  computation. 
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Fig.  7.12.  The  error  le | e( n)  of>  the  P_ver3ion  and  the  "pollution  free"  error. 

more  details  we  refer  to  Babuska,  Oh  [1986].  Similar  situation  occurs  also  for 
weaker  singularities  and  a  high  accuracy,  e.g.,  in  stresses.  The  accuracy  has  to  be 
achieved  in  practice  by  "filtering"  out  the  singularity  influence  by  using  few 
layers  (in  dependence  on  the  strength  of  the  singularity). 

8.  SOME  IMPLEMENT ATIONAL  ASPECTS 

The  h  and  h-p  versions  of  the  finite  element  method  give  a  large  flexibility  in 
the  change  of  degrees  of  the  elements  either  uniformly  or  element  by  element.  This 
has  to  be  respected  in  the  architecture  of  the  program.  One  of  the  essential  features 
is  the  use  of  the  hierarchical  shape  functions.  In  the  most  simple  setting  of  a  squar 
master  element  shown  in  Fig.  8.1.  The  hierarchical  shape  functions  are  as  follows: 

We  have  3  kinds  of  shape  functions 

0 

a)  The  modal  shape  functions  Ni(£,n). 


Fig.  8.1.  The  scheme  of  the  master  element. 


1 

b)  The  side  shape  functions  Nf(5,n). 

1  2 

c)  The  internal  shape  functions  N^Ct.n). 

The  index  above  describes  the  dimension  of  the  part  of  the  element  associated  with 
the  shape  functions. 

O 

a)  The  nodal  shape  functions  N^CC.n),  i  -  1 , . . . , H  are  the  usual  bilinear 

functions  associated  with  the  four  nodes. 

1 

b)  The  side  shape  functions  N^K.n),  i  *  1 . . . ,  4 ( p— 1 ) .  Let 


where  l  j ,  j  >  1  is  the  Legendre  polynomial.  The  the  side  shape  functions  asso¬ 
ciated  with  the  side  -1  <  C  <  1,  n  -  -1  are 


NjU.n)  -  P  (£)  ly1,  j  *  1 . p-1 

and  hence  we  have  4(p-1)  side  shape  functions. 

2  ? 

c)  The  internal  shape  functions  N^(£,n).  i  *  1 ..... ( p- 1 )  .  The  internal 
shape  functions  are 

NJikU.n)  -  Pj(OPk(n)  j  >  1.  k  >  i. 


The  number  of  the  internal  shape  functions  can  be  reduced  as  in  the  case  of  the 
serendipity  elements.  In  this  case  the  total  number  Q  of  all  of  the  shape  func¬ 
tions  is 

Q  -  4p  +  max(0,^(p-2) (p-3) 

as  used  in  the  program  PROBE.  See  Szabrf  [1985].  Probe  uses  1  i  p  $  8. 

The  shape  functions  are  called  hierarchical  because  the  increase  of  p  leads 
only  to  the  addition  of  new  shape  functions.  The  use  of  the  integrals  of  Legendre 


polynomials  has  essential  advantages  in  the  connection  with  the  round-off  error  con¬ 
trol  . 

The  hierarchical  type  of  this  elements  leads  to  the  construction  of  the  local 
stiffness  matrix  by  bordering  technique  when  p  increases.  The  same  structure  can 
be  achieved  for  the  global  stiffness  matrix.  This  character  then  allows  the  flex¬ 
ibility  to  change  the  degree  of  elements  and  its  use  for  the  quality  control  and 
adaptive  approaches. 

So  far  we  mentioned  only  the  square  elements.  Similar  hierarchical  structure 
can  also  be  used  for  the  triangular  elements. 

The  p-version  (and  h-p  version)  utilizes  large  elements.  Hence  it  is  important 
to  deal  properly  with  the  curved  boundaries.  It  seems  to  be  advantageous  (as  imple¬ 
mented  in  PROBE)  that  the  boundary  is  described  exactly  and  not  to  be  approximated, 
e.g.  by  polynomials.  This  is  achieved  by  using  blending  function  mappings.  This 
approach,  of  course,  leads  to  the  loss  of  the  important  rigid  body  motion  property 
of  the  isoparametric  elements.  Nevertheless,  if  the  boundaries  are  smooth,  e.g. 
piecewise  analytic,  then  for  higher  p  (say  p  -  3-4),  the  effect  of  this  loss  is  in 
practice  negligible  as  shown  theoretically  and  by  experience  with  PROBE.  In  addition, 
the  quality  control  indicates  also  this  possible  effect  (see  also,  Section  10). 

The  architecture  of  the  p-version  finite  element  program  has  essential  differ¬ 
ences  in  comparison  to  the  h-version.  We  mention  here  the  type  of  meshes  which  are 
used  which  influence  the  mesh  generator  and  many  other  aspects.  We  will  not  elaborate 
more . 

9.  THE  ASSESSMENT  OF  THE  COMPLEXITY  AND  THE  COMPUTATIONAL  COST 

The  cost  of  using  any  finite  element  code  consists  of  the  computer  cost  and  the 
human  cost.  The  relative  part  of  the  human  cost  growths  due  to  the  progress  in  the 
hardware  development  and  is  often  the  major  part  of  the  computational  analysis. 
Hence,  the  assessments  of  the  complexity  and  the  total  cost  of  the  method  is  very 
subtle  topic.  Some  of  the  aspects  will  be  also  discussed  in  Section  12  where  the  use 
of  the  h-p  version  and  the  experience  with  the  code  PROBE  in  the  industry  will  be 
addressed . 

The  computer  cost  of  running  the  code  consists  of  the  cost  of  the  basic  compu¬ 
tation  itself,  the  cost  of  the  assessment  of  the  quality  of  the  results  and  the  cost 
of  various  "bells"  and  "whistles"  and  the  "niceties"  for  the  user  which  any  commercial 
code  has  to  have.  The  computational  cost  depends  also  on  the  robustness  of  mesh 
design  criterion  and  selection  of  other  parameters  leading  to  the  desired  accuracy, 
and  how  effectively  the  user  is  able  to  make  the  proper  choice  of  these  parameters. 

The  computer  cost  is  the  cost  of  numerical  operations  and  the  10  cost  which  often 
can  be  the  main  cost.  This  makes  any  comparison  a  difficult  task.  It  seems  that 
one  of  the  best  ways  is  to  compare  various  codes  in  an  industrial  environment  (see 


In  general  we  have  to  relate  the  cost  of  the  achieved  (a-priori  given  by  the 
user)  accuracy  and  the  assessment  of  the  reliability  of  the  computed  data  (see  also 
Sections  10  and  11).  In  previous  sections  we  measured  the  effectivity  of  the  method 
by  the  relation  between  the  accuracy  measured  in  the  energy  norm  and  the  number  of 
degrees  of  freedom.  Of  course,  the  number  of  degrees  of  freedom  does  not  reflect 
some  aspects  of  computational  effort  as  the  computation  of  the  local  stiffness  ma¬ 
trices,  the  influence  of  the  sparsity  of  the  global  stiffness  matrix  and  hence  more 
detailed  analysis  both  theoretical  and  experimental  is  in  place.  The  only  count  of 
arithmetical  operations,  although  important,  is  not  sufficiently  realistic  basis  for 
the  computational  cost.  Hence,  very  likely  the  most  reliable  way  is  to  write  a  pro¬ 
gram  which  implements  the  method  and  analyze  its  performance.  To  this  end,  we  used 
part  of  the  program  PROBE  adjusted  to  given  purpose.  Although  the  results  depend  on 
the  programming,  the  computer,  compiler,  etc.,  the  assessment  based  on  such  an  anal¬ 
ysis  is  the  most  reliable  technique. 

The  major  parts  of  the  finite  element  computations  are 

a)  input,  treatment  of  the  geometry,  controls  of  the  input,  etc. 

b)  Computation  of  the  local  stiffness  matrices  and  load  vectors. 

c)  Assembly  and  elimination  procedure. 

d)  Postprocessing  inclusive  accuracy  assessment. 

We  will  address  here  only  parts  b)  and  c).  To  analyze  this  question,  let  us 
consider  the  mesh  shown  in  Fig.  9.1  with  r  $  m. 

We  mention  that  his  mesh  is  topologically  very  close  to  the  mesh  shown  in  Fig. 
7.6  (r  -  min(6,n+1),  m  -  max(6,n+1 ))  or  in  Fig.  7 .11.  The  cost  of  the  parts  b) 
and  c)  of  the  computations  dealing  with  the  examples  treated  in  Section  7  (Figs. 
7.5,  7.11)  is  essentially  the  same  as  when  dealing  with  the  mesh  of  Fig.  9.1  for 
properly  chosen  r  and  m. 


Fig.  9.1.  The  mesh  for  the  complexity  analysis. 

Let  us  now  count  the  number  of  operations  for  the  two  major  parts  b)  c)  men¬ 
tioned  above  for  the  elasticity  problem  (2  unknown  functions) 

a)  Computation  of  the  local  stiffness  matrix.  It  consists  of 

1)  Computation  of  the  values  of  the  shape  function  and  the  Jacobian  in 
one  Gaussian  point. 

2)  Computation  of  the  values  in  »  p2  Gaussian  points  and  the  sum- 


mation  for  obtaining  one  of  the  entries  of  the  local  stiffness  matrix. 

3)  Computations  of  all  entries  »  p^  of  the  local  stiffness  matrix. 

*0  Computations  of  mn  local  stiffness  matrix. 

7 

The  crude  operation  count  indicates  that  we  can  expect  s  «  mn  p  operations. 
Nevertheless,  this  count  is  deceptive  as  will  be  seen  below. 

The  number  of  the  shape  functions  is 

Q  -  [4p  +  max(0,^<p-2)(p-3))]2  (9.1) 

(the  serendipity  type  of  elements  are  used  and  hence,  for  IS  p  1  8  the 

second  term  is  not  too  large  and  the  number  of  the  shape  functions  is  closer 

p  D 

to  =  p  than  p  .  The  number  of  Gaussian  points  used  is  2[^]  +  2  where  [a] 

h 

denotes  the  integral  part  of  a.  A  careful  programming  leads  then  to  =  p 
operations  in  the  given  range.  In  Fig.  9.2  we  show  the  time  for  computations 
of  25  local  stiffness  matrices  (m  -  r  -  5)  in  dependence  on  p  in  the  dou¬ 
ble  logarithmic  scale  (time  units  on  VAX).  The  straight  line  with  the  shape 
4  confirms  the  rate  p1*  we  expected.  The  minor  "wabling"  of  the  time  curve 
is  caused  by  using  only  even  number  of  the  Gaussian  points, 
b)  The  assembly  and  elimination  form  m  >  r.  The  frontal  solver  has  been  used 
which  joints  the  elimination  and  the  assembly.  The  front  is  of  order  rp  (when 
properly  counting  (condensation  approach)  the  dealing  with  internal  shape  functions. 
Hence,  one  can  expect  the  time  is  of  order  «  (rp)2  mrp  =  mr^p^.  Fig.  9.3  shows  the 


Fig.  9.2.  The  computer  time  for  computation  of  25  local  stiffness  matrices. 


time  in  dependence  on  p  for  various  values  m  and  r.  The  slope  pJ  is  depicted 
in  the  figure  which  agree  with  the  expectation. 

Fig.  9.4  shows  the  time  in  dependence  on  m  with  m  =  r  and  different  p. 

The  slope  of  m^  is  shown  in  the  figure.  We  see  that  in  the  given  range,  the  rate 

o  li 

i3  closer  to  mJ  than  to  m  as  one  could  expect.  For  small  p  and  m  the  tests 
are  outside  of  the  asymptotic  range  which  had  to  be  expected. 

Analyzing  Figs.  9. 2-9. 4  and  the  experience,  we  see  that  the  number  of  opera¬ 
tions  for  the  h-p  version  with  number  of  layers  roughly  proportional  to  p  (see, 
e.g.  Fig.  7.9)  leads  to  number  of  operation  of  order  Na  where  1.5  £  a  £  2.5, 
where  with  p  increasing  a  increases.  The  conventional  h-ver3ion  with  p  « 

1,2  will  practically  never  lead  to  the  accuracy  of  order  of  1J  for  any  reasonable 
computational  effort  if  the  problem  is  not  very  simple.  Nevertheless,  for  small 
accuracy,  a  small  p  is  desirable  (see  Fig.  7.9).  Hence,  the  proper  design  of  the 
mesh  and  the  degree  of  elements  is  very  important.  This  will  be  discussed  in  the 
next  two  sections. 

The  method  is  very  robust  with  respect  to  the  mesh  design  when  some  basic  rules 
are  observed.  The  accuracy  can  be  obtained  by  using  higher  degree  without  changing 
the  mesh. 


Fig.  9.3.  The  time  of  the  assembly  and  the  elimination  in  dependence  on  p,  m,  r. 
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Fig.  9. *4.  The  time  of  the  assembly  and  the  elimination  for  r  -  m. 

10.  A-POSTERIORI  ASSESSMENT  OF  THE  QUALITY  OF  COMPUTED  DATA 

An  essential  part  of  the  computational  analysis  is  the  a-posteriorl  assessment 
of  the  quality  oi  the  obtained  results  of  interest.  (For  the  survey  of  the  main 
ideas  for  the  quality  assessment,  we  refer  to  Noor,  Babuska  [1986]).  Although  the 
energy  norm  is  basic  measure  of  accuracy,  it's  mostly  not  the  most  relevant  measure 
for  engineering  purposes. 

The  basic  and  practically  most  efficient  method  for  assessment  of  the  accuracy 
of  the  data  of  interest  is  comparing  computed  data  obtained  for  various  p  and  use 
of  the  extrapolation  method.  The  computed  data  can  be  of  various  type  inclusive 
various  equilibration  checks,  etc.  Although  PROBE  provides  many  reliability  checks, 
theoretically  the  a-posteriorl  quality  assessment  is  the  least  mathematically  under¬ 
stood  part  of  the  computation. 

The  assessment  of  the  error  measure  in  the  energy  form  is  the  basic  one  and 

p 

relatively  easiest  to  obtain  because  |e|E  -Eg  "Egg  where  Eg  (respectively  EF£) 
is  the  strain  energy  for  the  exact  (respectively,  finite  element)  solution  and  Efe 
increases  with  increase  of  p.  As  we  have  seen  earlier,  the  rate  of  convergence  of 
the  p-version  with  properly  designed  mesh  is  higher  than  the  algebraic  one.  Hence, 
one  can,  for  example  (as  implemented  in  PROBE)  assume  that 
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with  C  and  8  to  be  computed  from  3  consecutive  degrees  p  -  2,  p-1,p.  Because 

of  the  reverse  S  property  of  the  accuracy  mentioned  earlier  8  growths  with  p 
and  then  decreases  to  the  limiting  value  when  algebraic  rate  is  achieved.  The  de¬ 
sired  accuracy  should  be  practically  achieved  in  the  range  of  p  when  8  is  about 
maximal.  This  gives  an  indication  of  the  quality  of  the  mesh  for  given  accuracy. 

We  will  illustrate  the  approach  on  our  L-shape  problem  with  the  mesh  of  two 
layers  (n  «  2)  discussed  in  Section  7.  Table  10.1  gives  the  basic  data.  Column  5 
is  the  value  of  the  energy  norm  in  t  of  the  true  error.  Column  7  shows  the  relative 
error  in  %  (C.P.R.E)  when  computed  from  the  current  values  p  -  2,  p  -  1,p  of  the 
degrees.  Column  8  shows  the  error  (PRE)  computed  from  the  predicted  strain  energy 
for  p  «  6,7,8. 

Coming  back  to  Fig.  7.8  we  see  that  the  value  p  =  U  +■  5  which  is  indicated 
by  8  as  the  optimal  one  is  close  to  the  envelope  of  the  curves.  The  predicted 
error  CPRE  is  more  reliable  in  the  phase  when  8  is  increasings  than  decreasing 
which  is  obvious  because  of  (10.1).  Nevertheless,  as  a  whole  CPRE  is  very  reliable. 

Obviously,  PRE  is  very  accurate  for  lower  p. 

The  reliability  of  the  estimation  we  have  shown  above  is  closely  related  to  the 
monotone  behavior  of  the  error.  Nevertheless,  it  is  necessary  to  estimate  also 
other  values  of  interest.  We  shall  show  the  computation  of  the  stress  intensity 
factor  K  for  the  case  of  the  L-shaped  domain  analyzed  in  Section  7.  The  method  of 
computation  of  this  stress  intensity  was  developed  in  BabuSka,  Miller  [1984]  and 
Szabo,  BabuSka  [1986],  It  has  been  shown  that  this  method  gives  accuracy  which  is 
roughly  as  the  accuracy  of  the  strain  energy  (i.e.  the  square  of  the  energy  norm 
error).  Table  10.2  shows  the  values  of  the  stress  Intensity  factor  (normalized  to 
value  1).  For  more,  see  Szabo,  Babuska  [1986].  In  contrast  to  the  energy,  we  see 
that  the  error  changes  sign  and  Is  not  monotonic.  This,  of  course,  complicates  the 
extrapolation  procedure. 

Finally,  let  us  show  a  more  complicated  example  (where  the  exact  solution  is 
not  known)  of  the  cracked  panel  as  shown  in  Fig.  10.1  with  the  used  mesh. 

Table  10.3  shows  the  values  of  the  first  and  second  stress  intensity  factor 
Kj  and  Kjj. 

The  program  PROBE  has  built  in  additional  quality  tests  based  on  the  equilib¬ 
rium  testing  so  that  the  user  can  place  maximal  confidence  in  the  data  of  interest. 

Let  us  underline  that  because  of  the  hierarchic  structure  of  the  program  the 
computation  for  different  degrees  can  be  arranged  in  a  computationally  very  effec¬ 
tive  way. 

As  we  have  seen,  it  is  important  to  select  effectively  the  parameters  of  the  h- 
p  version,  i.e.  the  mesh  and  the  degree  of  the  elements  which  lead  to  the  desired 
accuracy  of  the  computed  data  of  interest.  This  can  be  done  by  a  feedback  (adap¬ 
tive)  approach  or  by  help  of  an  expert  system  or  by  the  combination  of  both. 
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TABLE  10.1 


3 

4 

5 

6 

7 

8 

le  1 E 

FE 

E’  FE 

1^* 

a 

C.P. R. E. 

P.R.E. 

3.88608797 

2.684  (-1) 

25.41 

_ 

_ 

25.41 

4.12483255 

2.971  (-2) 

8.45 

— 

— 

8.45 

4.14811501 

6.429  (-3) 

3-93 

1.93 

5.34 

3-91 

4.15265042 

1.893  (-3) 

2.13 

2.76 

2.02 

2.09 

4.15362580 

9.084  (-4) 

1.47 

3.05 

1 .01 

1 .41 

4.15297464 

5.695  (-4) 

1.37 

2.45 

.80 

1.10 

4.15413900 

4.052  (-4) 

.98 

1 .83 

.75 

.89 

4.15423777 

3.064  (-4) 

.85 

1.39 

.74 

.74 

4.15454422 

i 

i 


» 
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TABLE  10.2 


1 

2 

3 

P 

N 

K 

1 

41 

.  95268 

2 

119 

1.02177 

3 

209 

1.00250 

4 

335 

1 .00073 

5 

497 

0.99991 

6 

695 

0.99985 

7 

929 

0.99987 

8 

1199 

0.99990 

m 

OB 

1 .00000 
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TABLE  10.3 


1 

2 

3 

4 

p 

N 

KI 

KII 

1 

43 

.42259 

-.37480 

2 

125 

.55588 

-.25578 

3 

221 

.56161 

-.28951 

4 

355 

.59232 

-.28319 

5 

527 

.59825 

-.29398 

6 

737 

.60043 

-.28897 

7 

985 

.60119 

-.29196 

8 

1271 

.60132 

-.29042 

1 1 .  THE  EXPERT  SYSTEM  AND  THE  FEEDBACK  APPROACH 

Remark  11.1.  In  Gul ,  Babu3ka  [1986]  (Part  III)  we  developed  and  theoretically 
studied  a  feedback  system  which  leads  to  the  simultaneous  selection  of  the  mesh  and 
(nonuniform)  degrees  of  elements.  This  study  has  shown  that  the  adaptive  approach 
for  the  h-p  version,  based  on  similar  principles  as  for  the  h  version,  is  very  com¬ 
plicated  and  other  approach  has  to  be  likely  considered.  A 

The  feedback  system  for  the  p-version  with  uniform  p  (i.e.  for  fixed  me3h) 
can  be  relatively  easy  implemented.  We  have  seen  in  Sections  8  and  9  that  the  hier¬ 
archic  structure  of  the  p-version  allows  effectively  change  the  degree.  Hence,  the 
feedback  of  the  p-version  for  the  given  mesh  consists  of  successive  computation  for 
increasing  p  until  the  value  of  the  interest  will  have  admissible  accuracy.  The 
main  question  i3  what  is  the  ratio  of  the  cost  of  the  entire  computation  to  the  cost 
of  the  final  one.  Assuming  that  most  ineffective  computation  mode,  i.e.  independent 
computation  for  every  p  starting  from  p  -  1  and  that  the  cost  W(p)  of  the  com¬ 
putation  for  degree  p  is  W(p)  -  Cpa  the  ratio  for  (1  £  p  £  8)  for  a  -  3 
(respectively,  a  «  4)  ranges  from  1  to  2.5  (respectively,  2.14).  Hence,  a 
more  effective  computation  arrangement  exploiting  the  hierarchic  structure  leads  to 
the  ration  in  the  range  1.5-2  at  the  most.  This  ratio  is  very  comparable  with  the 
ratio  in  the  adaptive  approaches  based  on  the  h-version  for  the  energy  norm  accura¬ 
cy.  The  feedback  p-version  can  be  easily  used  for  any  accuracy  definition. 

Above  we  have  assumed  that  all  computations  starting  from  p  -  1  are  made  and 
that  the  mesh  is  given.  The  aim  is  obviously  to  construct  the  mesh  and  p  -  p0  so 
that  the  desired  accuracy  is  achieved  and  to  compute  only  the  solutions  for  p  £ 

Pq  (which  is  very  cheap  because  of  the  hierarchical  structure)  to  get  the  error 
estimates  based  on  the  extrapolation  procedure. 

Recently,  a  lot  of  progress  has  been  made  in  the  general  developing  of  the 
knowledge  based  systems  (expert  systems)  in  general  (see,  e.g.  Waterman  [1986], 
Hayes-Rogh,  Waterman,  Lenat  [1984])  and  various  attempt  has  been  made  to  apply  the 


system  for  computational  mechanics  problems  (see  e.g.,  Taig  [1986],  Fenves  [[1986]). 
In  Babuska,  Rank  [1986]  and  Rank,  Babuska  [1986]  and  expert  system  was  developed 
which  advises  the  user  how  to  use  the  h-p  version  in  an  effective  way.  One  of  the 
main  problems  is  to  develop  the  set  of  rules  which  is  then  utilized  in  the  system. 

The  finite  element  analysis  can  be  essentially  separated  into  two  parts,  the 
decision  phase  and  the  execution  phase,  which  carries  out  the  decision  although 
these  two  phases  are  often  interwoven.  The  adaptive  approach  combines  these  two 
phases  very  tightly  when  user  supplies  only  the  basic  data  and  the  program  takes 
care  about  the  entire  process.  The  approach  in  the  expert  system  mentioned  above, 
clearly  separates  the  phases  of  decision  and  execution  and  adds  the  part  of  the  ad¬ 
vise  and  communication.  It  is  using  theoretical  mathematical  results  mentioned  in 
previous  sections  in  combination  with  heuristics.  The  system  was  developed  for 
treating  polygonal  domains  and  accuracy  measured  in  the  energy  norm.  We  will  show 
now  an  example.  For  more  details,  we  refer  to  Babuska,  Rank  [1986]. 

Let  the  problem  be  the  plane  strain  of  the  "cracked  wrei.ch"  (symmetric)  domain 
shown  in  Fig.  11.1.  The  first  task  is  the  construction  of  the  geometry  mesh.  Fig. 
11.2  shows  it.  From  the  theory  of  the  h-p  version,  it  is  known,  that  areas  in  the 
neighborhood  of  the  corners  (entrant)  are  critical  and  have  to  be  jealt  with.  The 
user  is  advised  about  it  and  a  basic  mesh  which  separate  the  critical  areas  is  con¬ 
structed.  Fig.  11.3  shows  this  basic  mesh.  Program  then  computes  rough  solution  on 
basic  mesh,  computes  the  few  stress  Intensity  factors  for  every  corner  and  computes 
an  optimization  problem  leading  to  the  highst  accuracy  for  the  minimal  cost,  and 
gives  for  given  p  the  (optimal)  distribution  of  layers.  Based  on  this  computation 
the  expert  presents  the  user  with  the  relation  between  predicted  error  and  cost  ( N) 
for  the  optimal  mesh  computation.  Fig.  11. <1  shows  this  graph.  The  user  selects  the 
accuracy  and  the  program  constructs  the  mesh  for  final  solution.  Fig.  11.5  shows 
thi3  mesh  (not  all  layers  are  shown).  For  example,  the  mesh  for  the  error  1.8%  (re¬ 
spectively  1.3)  leads  to  the  predicted  optimal  degree  p  -  4 .  N  -  2101  (respec¬ 
tively  2517)  and  number  of  layers  shown  in  Table  11.1.  We  denoted  by  SL  (column  1) 
the  number  of  the  critical  area  as  shown  in  Fig.  11.1  and  by  L  the  suggested  number 
of  layers  in  this  area.  The  computation  on  this  mesh  leads  to  the  errors  3-1%  and 
2.6%,  respectively. 

We  see  that  the  predicted  errors  are  good  quality  although  not  perfect  what 
could  be  expected.  Nevertheless,  the  optimal  distribution  of  the  layers  was 
predicted  correctly. 

If  the  error  3.1%  is  not  acceptable,  then  p  can  be  increased  in  the  adaptive 
mode.  In  our  case,  this  would  mean  to  increase  N  to  3059  which  will  lead  to  the 
error  1.9%.  (The  other  possibility  is  to  construct  the  optimal  mesh  for  higher 
accuracy . ) 

Let  us  mentioned  that  the  expert  system  we  outlined  cost  about  10-20%  of  the 


entire  computer  cost. 
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Fig.  11.6  shows  the  actual  error  for  the  meshes  and  degrees  shown  in  Fig.  11.4. 
Finally,  let  us  underline  that  the  expert  system  and  the  problem  of  the  optimal 
mesh  design  is  in  a  very  beginning  stage. 
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Fig.  11.6.  Actual  error  for  cracked  wrench  problem. 

12.  INDUSTRIAL  EXPERIENCE  WITH  THE  p  AND  h-p  VERSIONS  OF  THE  FINITE  ELEMENT  METHOD 

As  was  said  in  Section  1  the  p-version  of  the  finite  element  method  was  first 
introduced  in  USA  when  the  program  FIESTA  was  released  -  Peano,  Walker  [1981].  The 
first  release  of  FIESTA  had  elastostatic  capabilities  for  solving  three  dimensional 
problems.  The  p-extension  (i.e.,  Increase  of  p)  was  possible  through  the  specif i 
cation  six  group  of  hierarchic  basis  function  corresponding  to  polynomial  degrees 
ranging  from  1  to  4.  The  program  had  appealed  to  engineers  mainly  because  it  al¬ 
lowed  the  use  of  coarse  meshes  and  therefore  simplified  data  preparation  tasks. 
Visual  checking  of  the  mesh  in  three  dimensions  is  generally  very  difficult  and 
often  impossible  when  the  h-version  is  used. 

The  industrial  users  are  usually  far  more  concerned  about  human  time  than  ma¬ 
chine  cycle  requirements  but  thus  greatest  concern  is  the  elapsed  real  time  between 
problem  definition  and  delivery  of  a  reliable  solution.  This  is  because  the  costs 
associated  with  numerical  analysis  are  usually  negligible  as  compared  with  the  con¬ 
sequences  of  not  being  able  to  deliver  correct  solutions  on  time  within  the  project 
schedule  prescribed  by  management. 

The  first  release  of  the  program  PROBE  was  introduced  to  the  aerospace  industry 
in  June  1985  by  Noetic  Technologies  Corporation,  a  St.  Louis  firm.  This  release  had 
capabilities  in  linear  plane  elasticity  only.  The  second  release  also  has  capabil¬ 
ities  in  linear  two  dimensional  heat  transfer,  axisymmetric  and  planar  thermoelas- 
.  ticity.  The  first,  release  .was.of.  interest  mainly  to^airframe  designers,  whereas  .the 
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second  release  is  also  of  interest  of  engine  designers. 

The  program  PROBE  was  extensively  tested  and  is  used  very  successfully  in  the 
industry,  e.g.,  in  the  analysis  of  orthotropic  panels,  various  cut-cuts  fracture 
mechanics,  problems  involving  laminated  composites  where  the  peeling  stress  was  of 
primary  interest,  etc. 

Engineer  users  consider  various  features  to  be  of  important  for  them  and  judge 
the  programs  and  methods  accordingly.  The  features  are  as  follows: 

1)  Increased  level  of  confidence  in  the  computation.  The  usual  observation  of 
the  convergence  of  the  data  of  interest  (not  only  the  energy  norm)  and  various  addi¬ 
tional  tests  as  equilibrium  checks  (having  physical  meaning),  etc.  give  the  user  the 
confidence  in  the  computational  analysis.  (For  some  comments,  see  Section  8.) 

2)  Lower  human  time  requirement  due  to  simplicity  and  flexibility  of  input. 

In  an  industrial  test  30-40  fold  saving  in  human  time  were  reported  in  comparison 
with  conventional  finite  element  technology  (Barhart,  Eisemann  [1986]). 

3)  To  get  rapid  convergence  and  flexibility,  e.g.,  when  large  aspect  ration 
have  to  be  used  as  e.g.  in  the  case  of  composite  joints. 

4)  The  flexibility  of  mesh  design,  e.g.  in  fracture  mechanics  following  the 
crack  growth,  series  of  short  cracks,  optimal  design  problem,  etc.  (e.g.  Schiermein, 
Szabo  [1986]). 

5)  Easy  learning  and  robust  performance. 

The  p  and  h-p  version  is  very  well  responding  to  needs  of  this  type.  Never¬ 
theless,  further  new  features  as  dealing  with  plate  shells,  fabricated  plates,  e 
dimensional  problems,  nonlinear  analysis  need  to  be  developed. 

13.  RELATION  TO  SOME  OTHER  METHODS 

The  p  and  h-p  versions  of  the  finite  element  method  has  a  relation  to  todays 
form  of  the  spectral  method.  The  spectral  method  (see,  e.g.,  Gottlieb,  D.  D., 
Hussaini,  Orszag,  [1984])  expands  the  solution  of  the  problem  in  high-order  Fourier 
or  polynomial  series,  the  coefficients  of  which  are  determined  by  weighted-residual 
projections.  The  spectral  method  is  used  in  the  fluid  mechanics,  for  example,  in 
the  problem  of  transitional  and  low  Reynold  number,  turbulent  incompressible  fluid 
flow  in  simple  domains.  Recently,  the  research  of  the  spectral  method  is  focused  on 
extensions  to  more  complex  domains.  The  spectral  method  uses  Galerkin  type  approach 
and  is  assuming  that  the  conditions  (2.5)  holds.  The  spectral  method  is,  with  the 
use  of  polynomial  approximation,  closely  related  to  the  p-version  with  differences 
in  applications  and  implementations.  See  Patera  [1986]  and  the  references  given 
there.  The  pseudospectral  method  can  be  viewed  as  the  Galerkin  method  with  the  nu¬ 
merical  integration  technique. 

The  spectral  method  was  traditionally  analyzed  in  the  context  with  smooth  solu¬ 
tions.  Only  recently  the  needs  of  copping  with  the  singularities  and  geometries  in 
the  spirit  of  the  h-p  version  is  seen. 


The  p-version  is  also  related  to  the  Global  Element  Method  by  Delves  Phillips 
[1980]. 

The  p  and  h-p  versions  of  the  finite  element  method  is  relatively  well  devel¬ 
oped  and  taylored  for  the  needs  of  general  problems  in  structural  mechanics  as  we 
have  been  seen  in  previous  sections.  The  family  of  spectral  methods  which  is  ap¬ 
plied  in  the  field  of  fluid  mechanics  has  various  aspects  which  are  related  to  the  p 
and  h-p  versions  discussed  in  this  paper. 
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